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ABSTRACT 

This  report  forms  the  user’s  guide  for  Version  3.2  of  QPSOL,  a  set  of  Fortran  subroutines  designed 
to  locate  the  minimum  value  of  a  quadratic  function  subject  to  linear  constraints  and  simple  upper 
and  lower  bounds.  If  the  quadratic  function  is  convex,  a  global  minimum  »  found;  otherwise,  a 
local  minimum  is  found.  The  method  used  is  most  efficient  when  many  constraints  or  bounds  are 
active  at  the  solution.  QPSOL  treats  the  Hessian  and  general  constraints  as  dense  matrices,  and 
hence  is  not  intended  for  large  sparse  problems. 

This  document  replaces  the  previous  user’s  guide  of  July  1083. 


t  QPSOL  is  available  from  the  Office  of  Technology  Licensing,  105  Encina  Hall,  Stanford  University, 
Stanford,  California,  94305. 

The  material  contained  in  this  report  is  based  upon  research  supported  by  the  U.S.  Department  of 
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1.  PURPOSE 


QPSOL  is  a  collection  of  Fortran  subroutines  designed  to  solve  the  quadratic  programming 
(QP)  problem  —  the  minimization  of  a  quadratic  function  subject  to  a  set  of  linear  constraints 
on  the  variables.  The  problem  is  assumed  to  be  stated  in  the  following  form: 


QP 

minimize 

cTx  +  \xtHx 

xeR” 

2 

subject  to 

1  £  {£} £  “■ 

where  e  is  a  constant  n-vec tor  and  //  is  a  constant  n  X  n  symmetric  matrix;  note  that  II  is  the 
Hessian  matrix  (matrix  or  second  partial  derivatives)  of  the  quadratic  objective  function.  The 
matrix  A  is  m  X  n,  where  m  may  be  zero;  A  is  treated  as  a  dense  matrix. 

The  constraints  involving  A  will  be  called  the  general  constraints.  Note  that  upper  and  lower 
bounds  are  specified  for  all  the  variables  and  for  all  the  general  constraints.  The  form  of  QP 
allows  Tull  generality  in  specifying  other  types  of  constraints.  In  particular,  an  equality  constraint 
is  specified  by  setting  li  —  If  certain  bounds  arc  not  present,  the  associated  elements  of  l  or 
u  can  be  set  to  special  values  that  will  be  treated  as  — oo  or  +oo. 

The  user  must  supply  an  initial  estimate  of  the  solution  to  QP,  and  a  subroutine  that  computes 
the  product  Ux  for  any  given  vector  x.  Some  typical  examples  of  this  subroutine  arc  included 
with  the  QPSOL  package.  There  is  no  restriction  on  II  apart  from  symmetry.  If  If  is  positive 
definite  or  positive  semi-definite,  QPSOL  will  obtain  a  global  minimum;  otherwise,  the  solution 
obtained  will  be  a  local  minimum  (which  may  or  may  not  be  a  global  minimum).  If  II  is  defined 
as  the  zero  matrix,  QPSOL  will  solve  the  resulting  linear  programming  (LP)  problem;  however, 
this  can  be  accomplished  more  ellicicntly  by  setting  a  logical  variable  in  the  call  of  subroutine 
QPSOL  (see  the  parameter  LP  in  Section  4),  or  by  using  the  LPSOL  package. 

QPSOL  allows  the  user  to  provide  the  indices  of  the  constraints  that  arc  believed  to  be 
satisfied  exactly  at  the  solution.  This  facility,  known  as  a  warm  start,  can  lead  to  significant 
savings  in  computational  olTort  when  solving  a  sequence  of  related  problems.  For  example,  the 
NPSOL  package  or  Gill  el  al.  (1984b)  uses  this  feature  in  a  sequential  quadratic  programming 
method  for  nonlincarly  constrained  optimization. 

The  quantity  of  output  is  controlled  by  the  user  (see  the  parameter  MSGLVL  discussed  in 
Section  4).  The  QPSOL  package  contains  approximately  6000  lines  of  ANSI  (1966)  Standard 
Fortran,  of  which  44%  arc  comments. 
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2.  DESCRIPTION 

The  method  used  to  solve  QP  is  an  active-set  null-space  method,  and  is  described  in  detail 
in  Gill  ct  al.  (19840);  a  closely  related  method  is  given  in  Gill  and  Murray  (1978).  The  main 
features  of  the  method  are  presented  here.  Where  possible,  explicit  reference  is  made  to  the  names 
of  variables  that  arc  parameters  of  subroutine  (JPSOL  or  are  mentioned  in  the  printed  output. 

The  method  has  two  distinct  phases.  In  the  first  (the  LP  phase),  an  iterative  procedure  is 
carried  out  to  determine  a  feasible  point.  In  this  context,  feasibility  is  defined  by  a  user-provided 
array  FEATOL;  the  j- th  constraint  is  considered  satisfied  if  its  violation  docs  not  exceed  FEATOL(j) 
(see  the  discussion  of  FEATOL  in  Section  4.)  The  second  phase  (the  QP  phase)  generates  a  sequence 
of  feasible  iterates  in  order  to  minimize  the  quadratic  objective  function.  In  both  phases,  a  subset 
of  the  constraints  called  the  working  set  —  is  used  to  define  the  search  direction  at  each 
iteration;  typically,  the  working  set  includes  constraints  that  are  satisfied  “exactly”  (to  within 
the  corresponding  tolerances  in  the  FEATOL  array). 

We  now  briefly  describe  a  typical  iteration  in  the  QP  phase.  Let  xk  denote  the  estimate  of 
the  solution  at  the  &-th  iteration;  the  next  iterate  is  defined  by 

Xk+l  =  Xk  +  OtkPk, 

where  pk  is  an  n-dimcnsional  search  direction  and  ak  is  a  scalar  step  length.  Assume  that 
the  working  set  contains  tk  linearly  independent  constraints,  and  let  Ck  denote  the  matrix  of 
coefficients  of  the  bounds  and  general  constraints  in  the  current  working  set. 

(.cl  Zk  denote  a  matrix  whose  columns  form  a  basis  Tor  the  null  space  of  Ck,  so  that  CkZk  — 
0.  (Note  that  Zk  has  nz  columns,  where  nx  =  n  —  tk.)  The  vector  Zk(c  +  //xt)  is  called  the 
projected  gradient  at  x*.  If  the  projected  gradient  is  zero  at  x*  (i.c.,  x*  is  a  constraint!  stationary 
point  in  the  subspacc  defined  by  Zk),  Lagrange  multipliers  X*  are  defined  as  the  solution  of  the 
compatible  overdetermined  system 

CTk\k=c  +  Hxk .  (1) 

The  Lagrange  multiplier  X  corresponding  to  an  inequality  constraint  in  the  working  set  is  said  to 
be  optimal  if  X  <  0  when  the  associated  constraint  is  at  its  upper  bound,  or  iT  X  >  0  when  the 
associated  constraint  is  al  its  lower  bound.  If  a  multiplier  is  non-optitnal,  the  objective  function 
can  be  reduced  by  deleting  the  corresponding  constraint  (with  index  KDEL)  from  the  working  set. 

If  the  projected  gradient  al  xk  is  nonzero,  the  search  direction  pk  is  defined  as 

Pk  =  ZkPw,  (2) 

where  pM  is  an  n*-vcctor.  In  effect,  the  constraints  in  the  working  set  arc  treated  as  equalities,  by 
constraining  p*  to  fie  within  the  subspace  of  vectors  orthogonal  to  the  rows  of  Ck.  This  definition 
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ensures  that  CkPk  —  0,  and  hence  the  values  of  the  constraints  in  the  working  set  are  not  altered 
by  any  move  along  p*. 

The  vector  pM  is  obtained  by  solving  the  equations 
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ZTkHZkpM  =  -ZTk{c  +  Hxk).  (3) 

(The  matrix  ZkHZk  is  called  the  projected  Hessian  matrix.)  If  the  projected  Hessian  is  positive 
definite,  the  vector  defined  by  (2)  and  (3)  is  the  step  to  the  minimum  of  the  quadratic  function 
in  the  subspace  defined  by  Zk. 

If  the  projected  llcssian  is  positive  definite  and  Xk  +  pt  is  feasible,  otfc  will  be  taken  as  unity. 
In  this  case,  the  projected  gradient  at  xk+l  will  be  icro  (sec  the  variable  NORM  ZTG  in  the  output 
from  QPSOL),  and  Lagrange  mu lli pliers  can  be  computed  (bcc  (1)).  Otherwise,  a*  is  set  to  the 
step  to  the  “nearest"  constraint  (with  index  KADD),  which  is  added  to  the  working  set  at  the  next 
iteration. 

The  matrix  Zk  is  obtained  from  the  TQ  factorisation  of  Ck,  in  which  Ck  is  represented  as 


CkQ  —  (  0  Tk),  (4) 

where  Tk  is  reverse-triangular  (sec  Gill  ct  a f.,  1984a).  It  follows  from  (4)  that  Zk  may  be  taken 
as  the  first  n,  columns  of  Q.  If  the  projected  Hessian  is  positive  definite,  (3)  is  solved  using  the 
Cholcsky  factorisation 

Z\HZk  =  /£«*, 

where  lik  is  upper  triangular.  These  factorisations  arc  update d  as  constraints  enter  or  leave  the 
working  set.  The  update  procedures  are  described  in  detail  in  Gill  ct  a I.  (1984a). 

An  important  feature  of  QPSOL  is  the  treatment  of  indefiniteness  in  the  projected  Hessian. 
If  the  projected  Hessian  is  positive  definite,  it  inay  become  indefinite  only  when  a  constraint  is 
deleted  from  the  working  set.  In  this  case,  a  temporary  modification  (or  magnitude  HESS  MOD)  is 
added  to  the  last  diagonal  clement  or  the  Cholcsky  factor.  Once  a  modification  has  occurred,  no 
further  constraints  are  deleted  from  the  working  set  until  enough  constraints  have  been  added  so 
that  the  projected  Hessian  is  again  positive  definite,  ir  problem  QP  has  a  finite  solution,  a  move 
along  the  direction  obtained  by  solving  (3)  with  the  modified  Cholesky  factor  must  encounter  a 
constraint  that  is  not  already  in  the  working  set. 

In  order  to  resolve  indcfinitcncss  in  this  way,  we  must  ensure  that  the  projected  Hessian  is 
positive  definite  at  the  first  iterate  in  the  QP  phase.  Given  the  nM  X  n,  projected  Hessian,  a 
step-wise  Cholcsky  factorisation  is  performed  with  symmetric  interchanges  (and  corresponding 
rearrangement  of  the  columns  of  Z),  terminating  if  the  next  step  would  cause  the  matrix  to 
become  indefinite.  This  determines  the  largest  possible  positive-definite  principal  submatrix  of 
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the  (permuted)  projected  Hessian.  If  n„  steps  of  the  Cholcsky  factorization  have  been  successfully 
completed,  the  relevant  projected  Hessian  is  an  n„  X  nM  positive-definite  matrix  ZTIf  X,  where 
X  comprises  the  first  nH  columns  of  X.  The  quadratic  function  will  subsequently  be  minimized 
within  subspaccs  of  reduced  dimension  until  the  full  projected  Hessian  is  positive  definite. 

Several  strategics  arc  used  to  control  ill-conditioning  in  the  working  set.  One  such  strategy 
is  associated  with  the  FEATOL  array.  Allowing  the  ji-th  constraint  to  be  violated  by  as  much  as 
FEATOL(y)  often  provides  a  choice  of  constraints  that  could  be  added  to  the  working  set.  When 
a  choice  exists,  the  decision  is  based  on  the  conditioning  of  the  working  set.  Negative  steps  arc 
occasionally  permitted,  since  i*  may  violate  the  constraint  to  be  added. 
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8.  SPECIFICATION 


SUBROUTINE  QPSOL(  ITMAX,  MSGLVL,  N, 

NCLIN,  NCTOTL ,  NROWA,  NROWH.  NCOLH, 

BIGBND ,  A,  BL.  BU,  CVEC.  FEATOL,  HESS,  QPHESS, 
COLD,  LP,  ORTHOG,  ISTATE,  X, 

INFORM,  ITER,  OBJ.  CLAMDA, 

IV.  LENIN,  V,  LENV  ) 


LOGICAL 

EXTERNAL 

INTEGER 

INTEGER 

REAL 

REAL 


COLD.  LP,  ORTHOG 
QPHESS 

ITMAX.  MSGLVL.  N.  NCLIN.  NCTOTL, 

NROWA.  NROWH.  NCOLH.  INFORM,  ITER.  LENIN,  LENV 
ISTATE (NCTOTL).  IW (LENIN) 

BIGBND.  OBJ 

A (NROWA, N) ,  BL (NCTOTL) ,  BU (NCTOTL) ,  CVEC(N), 
FEATOL (NCTOTL).  HESS (NROWH.  NCOLH).  X(N) , 
CLAMDA (NCTOTL) ,  W(LENW) 


Note:  Here  and  elsewhere,  the  specification  of  a  parameter  as  REAL  should  be  interpreted 
working  precision,  which  may  be  DOUBLE  PRECISION  in  some  circumstances. 
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4.  INPUT  PARAMETERS 


ITMAX  is  an  upper  bound  on  the  number  of  iterations  to  be  taken  during  the  LP  phase  or  the 
QP  phase. 

liSGLVL  indicates  the  amount  of  intermediate  output  desired.  The  printout  is  described  in 
Section  9.  All  output  is  written  to  the  file  number  NOUT  (sec  subroutine  UCHPAR  in 
Section  11).  For  MSGLVL  >  10,  each  value  of  liSGLVL  includes  the  output  from  all  lower 
values.  The  printout  corresponding  to  each  value  of  liSGLVL  is  defined  as  follows: 


MSGLVL 

Definition 

0 

No  output. 

1 

The  final  solution  only. 

5 

One  brief  line  of  output  for  each  constraint  addition  or  deletion  (no 
printout  of  the  final  solution). 

>  io 

The  final  solution  and  one  brief  line  of  output  for  each  constraint 
addition  or  deletion. 

>  15 

At  each  iteration,  X,  ISTATE,  and  the  indices  of  the  free  variables  (i.c., 
the  variables  not  currently  held  on  a  bound). 

>  20 

At  each  iteration,  the  Lagrange  multiplier  estimates  and  the  general 
constraint  values. 

>  30 

At  each  iteration,  the  diagonal  elements  of  the  matrix  T  associated 
with  the  TQ  factorisation  of  the  working  set,  and  the  diagonal  ele¬ 
ments  of  the  Cholcaky  factor  It  or  the  projected  Ilessian. 

>  80 

Debug  printout. 

99 

The  arrays  CVEC  and  HESS. 

N  is  the  number  of  variables  (i.c.,  the  dimension  ofX  ).  N  must  be  positive. 

NCLIN  is  the  number  or  general  linear  constraints  in  the  problem  (NCLIN  may  be  zero). 

NCTOTL  must  be  set  to  N  +  NCLIN. 

NROVA  is  the  declared  row  dimension  of  A  (NROWA  must  be  at  least  1  and  at  least  NCLIN). 

NROVH  is  the  declared  row  dimension  of  the  array  HESS  (NROWH  must  be  at  least  1). 

MCOLH  is  the  declared  column  dimension  of  the  array  HESS  (NCOLH  must  be  at  least  1). 
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BIGBND  is  a  positive  real  variable  whose  magnitude  denotes  an  “infinite"  component  of  t  and 
u.  Any  upper  bound  greater  than  or  equal  to  BIGBND  will  be  regarded  as  plus  infinity 
(and  similarly  for  a  lower  bound  less  than  or  equal  to  —BIGBND). 

A  is  a  real  array  of  declared  dimension  (NR0WA,N).  The  t-th  row  of  A  contains  the  coefficients 

of  the  t-th  general  constraint,  t  =  1  to  NCLIN.  If  NCLIN  is  zero,  A  is  not  accessed. 

BL  is  a  real  array  of  dimension  NCTOTL  that  contains  the  lower  bounds  for  all  the  constraints, 

in  the  following  order  (which  is  also  observed  for  BU,  ISTATE,  and  C LAUDA):  the  first 
N  elements  of  BL  contain  the  lower  bounds  on  the  variables;  if  NCLIN  >  0,  the  next 
NCLIN  elements  of  BL  contain  the  lower  bounds  for  the  general  linear  constraints.  In 
order  for  the  problem  specification  to  be  meaningful,  it  is  required  that  BL(j')  <  BU(j) 
for  all  j.  To  specify  a  non-existent  lower  bound  (i.e.,  tj  =  —  oo),  the  value  used  must 
satisfy  BL(j)  <  —BIGBND.  To  specify  the  j-th  constraint  as  an  equality,  the  user  must 
set  BL(j)  sb  BU(j')  =  j3,  say,  where  |/3|  <  BIGBND. 

BU  is  a  real  array  of  dimension  NCTOTL  that  contains  the  upper  bound  '  all  the  con¬ 

straints,  in  the  same  order  described  above  under  BL.  To  specify  a  no  .tent  upper 
bound  (i.e.,  uy  =  oo),  the  value  used  must  satisfy  BU(j)  >  BIGBND. 

CVEC  is  a  real  array  of  dimension  N  containing  the  coefficients  of  the  linear  term  of  the 
objective  function  (the  vector  c  in  problem  QP). 

FEATOL  is  a  real  array  of  dimension  NCTOTL  containing  positive  tolerances  that  define  the 
maximum  permissible  violation  in  each  constraint  in  order  for  a  point  to  be  considered 
feasible,  i.e.  constraint  j  is  considered  satisfied  ir  its  violation  does  not  exceed  FEATOL(y). 
Note  that  FEATOL(j)  is  a  bound  on  the  absolute  acceptable  violation.  For  example,  if  the 
data  defining  the  constraints  arc  of  order  unity  and  arc  correct  to  about  6  decimal  digits, 
it  would  be  appropriate  to  choose  FEATOL(j)  as  10-6  for  all  relevant  j.  In  genera),  the 
elements  of  FEATOL  should  be  chosen  as  the  largest  possible  acceptable  values,  since  the 
algorithm  of  QPSOL  becomes  less  likely  to  encounter  difficulties  with  ill-conditioning 
and  degeneracy  as  the  components  or  FEATOL  increase.  A  warning  message  is  printed 
if  any  component  or  FEATOL  is  less  than  machine  precision;  the  user  must  not  set  any 
component  or  FEATOL  to  zero.  A  detailed  discussion  or  FEATOL  is  given  in  Gill  ct  a/. 
(1984c). 

HESS  is  a  real  array  of  declared  dimension  (NROWH,  NCOLH)  that  may  be  used  to  store  the 
Ilessian  matrix  fl  of  problem  QP  if  desired.  The  elements  of  HESS  arc  accessed  only  by 
the  subroutine  QPHESS;  thus  HESS  is  not  accessed  if  LP  is  .TRUE.  In  some  cases,  the  user 
need  not  use  HESS  to  store  //  explicitly  (sec  the  specification  of  QPHESS  betow). 
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QPHGSS  is  the  name  of  a  subroutine  that  defines  the  Ilessian  matrix.  QPHESS  must  be  declared  as 
EXTERNAL  in  the  routine  that  calls  QPSOL.  QPHESS  is  not  accessed  if  the  logical  variable 
LP  is  .TRUE,  (see  the  description  below  of  LP).  The  user  has  considerable  flexibility  in 
coding  QPHESS  because  the  algorithm  of  Ql’SOL  requires  only  the  product  of  II  and  a 
vector;  the  elements  of  the  matrix  //  need  not  be  defined  explicitly.  Several  examples  of 
QPHESS  are  provided  in  order  to  demonstrate  some  of  the  alternatives.  The  specification 
of  QPHESS  is: 

SUBROUTINE  QPHESS (  N,  NROWH,  NCOLH,  JTHCOL,  HESS,  X,  HX  ) 


INTEGER 


N,  NROWH,  NCOLH.  JTHCOL 

HESS (NROWH,  NCOLH).  X(N) ,  HX(N) . 


The  actual  parameters  N,  NROWH,  NCOLH  and  HESS  input  to  QPHESS  will  always  be  the 
same  Fortran  variables  and  arrays  as  those  input  to  QPSOL.  They  must  no t  be  altered 
by  QPHESS. 

For  a  given  vector  x  (the  array  X),  the  array  HX  must  contain  the  product  Ilx  on 
exit  from  QPHESS. 

The  input  parameter  JTHCOL  is  included  to  allow  flexibility  Tor  the  user  in  the 
special  situation  when  *  is  the  j-th  coordinate  vector  (i.c.,  the  j-lh  column  of  the 
identity  matrix).  This  may  be  of  interest  because  the  product  Ilx  is  then  the  j-th 
column  of  II,  which  can  sometimes  be  computed  very  efficiently.  The  user  may  code 
QPHESS  to  take  advantage  of  this  case.  If  JTHCOL  =  j,  where  j  >  0,  HX  should  contain 
column  JTHCOL  of  //,  and  hence  special  code  may  be  included  in  QPHESS  to  test  JTHCOL 
if  desired.  However,  special  code  is  not  necessary,  since  the  vector  X  always  contains 
column  JTHCOL  of  the  identity  matrix  whenever  QPHESS  is  called  with  JTHCOL  >  0. 

In  some  cases,  it  may  be  desirable  to  use  a  one-dimensional  array  to  transmit  data 
or  workspace  to  QPHESS;  HESS  should  then  be  declared  as  REAL  HESS(NROWH),  and  the 
parameter  NCOLH  must  be  I.  (This  device  is  used  Tor  the  example  subroutines  QPHES4 
and  QPHES6  in  the  QPSOL  package,  to  economize  on  storage.) 

In  other  situations,  it  may  be  desirable  to  compute  // x  without  accessing  HESS — 
for  example,  if  II  is  sparse  or  has  special  structure.  (This  is  illustrated  in  the  subroutine 
QPHES1  in  the  QPSOL  package.)  The  parameters  HESS,  NROWH  and  NCOLH  may  then  refer 
to  any  convenient  array. 

When  MSGLVL  =  99,  the  (possibly  undefined)  contents  of  HESS  will  be  printed, 
except  if  NROWH  and  NCOLH  are  both  1.  Also  printed  arc  the  results  of  calling  QPHESS 
with  JTHCOL  =  1,2, _ ,  M. 

COLD  is  a  logical  variable  that  indicates  whether  the  user  wishes  to  specify  the  initial  working 
set.  In  general,  COLD  should  be  set  to  .TRUE,  for  the  first  call  of  QPSOL,  and  the 
initial  working  set  will  then  be  selected  by  QPSOL.  However,  if  a  good  estimate  of  the 
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LP 


ORTHOG 


initial  working  set  is  available  —  for  example,  when  QPSOL  is  called  repeatedly  to  solve 
related  problems  —  it  may  be  advantageous  to  set  COLD  to  .FALSE,  after  the  first 
call.  When  COLD  is  .FALSE.,  the  user  must  define  every  element  of  ISTATE  (see  the 
description  of  ISTATE  for  the  meaning  of  each  possible  value).  QPSOL  will  override  the 
user’s  specification  of  ISTATE  if  necessary,  so  that  a  poor  choice  of  the  working  set  will 
not  cause  a  fatal  error. 

is  a  logical  variable.  If  LP  is  .FALSE.,  QPSOL  will  solve  the  specified  quadratic  program- 
filing  problem.  If  LP  is  .TRUE.,  QPSOL  will  treat  H  as  zero  and  solve  the  resulting  linear 
program;  in  this  case,  parameters  HESS  and  QPHESS  will  not  be  accessed. 

is  a  logical  variable  that  indicates  whether  orthogonal  transformations  are  to  be  used 
in  computing  and  updating  the  TQ  factorization  of  the  working  set 

AQ  =  (o  r), 

where  A  is  a  submatrix  of  A  and  T  is  reverse-triangular.  If  ORTHOG  is  .TRUE.,  the  TQ 
factorization  is  computed  using  Householder  reflections  and  plane  rotations,  and  the 
matrix  Q  is  orthogonal.  IT  ORTHOG  is  .FALSE.,  stabilized  elementary  transformations  are 
used  to  maintain  the  factorization,  and  Q  is  not  orthogonal.  A  rule  of  thumb  in  making 
the  choice  is  that  orthogonal  transformations  require  more  work,  but  provide  greater 
numerical  stability.  Thus,  we  recommend  setting  ORTHOG  to  .TRUE.  iT  the  problem 
is  reasonably  small  or  the  active  set  is  ill-conditioned.  Otherwise,  setting  ORTHOG  to 
.FALSE,  will  often  lead  to  a  reduction  in  solution  time  with  negligible  loss  of  reliability. 
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ISTATE  is  an  integer  array  of  dim  •’sion  NCTOTL  that  indicates  the  status  of  every  constraint 
with  respect  to  the  working  set.  The  ordering  of  ISTATE  is  the  same  as  that  described 
above  under  BL,  i.e.,  the  first  N  components  of  ISTATE  refer  to  the  upper  and  lower 
bounds  on  the  variables,  and  components  N  +  1  through  N  +  NCLIN  refer  to  the  upper 
and  lower  bounds  on  Ax.  The  significance  of  each  possible  value  of  ISTATE(j)  is  as 
follows: 


ISTATE(j) 

-2 

-1 

0 

1 


2 


3 


Meaning 

The  constraint  violates  its  lower  bound  by  more  than  FEATOL(j).  This 
value  of  ISTATE  cannot  occur  after  a  feasible  point  has  been  found. 
The  constraint  violates  its  upper  bound  by  more  than  FEATOL(ji).  This 
value  of  ISTATE  cannot  occur  after  a  feasible  point  has  been  found. 
The  constraint  is  not  in  the  working  set.  Usually,  this  means  that  the 
constraint  lies  strictly  between  its  bounds. 

This  inequality  constraint  is  included  in  the  working  set  at  its  lower 
bound.  The  value  of  the  constraint  is  within  FEATOL(j)  of  its  lower 
bound. 

This  inequality  constraint  is  included  in  the  working  set  at  its  upper 
bound.  The  value  of  the  constraint  is  within  FEATOL(y)  or  its  upper 
bound. 

The  constraint  is  included  in  the  working  set  as  an  equality.  This  value 
or  ISTATE  can  occur  only  when  BL(j)  =  BU(j).  The  corresponding 
constraint  is  within  FEATOL(j)  of  its  required  value. 


If  COLD  =  .TRUE.,  ISTATE  need  not  be  set  by  the  user.  However,  when  COLD  is 
.FALSE.,  every  element  of  ISTATE  must  be  set  to  one  of  the  values  given  above  to  define 
a  suggested  initial  working  set  (which  will  be  changed  by  QPSOL  if  necessary).  The  most 
likely  values  are: 

ISTATE(j)  Meaning 


0  The  corresponding  constraint  should  not  be  in  the  initial  working  set. 

1  The  constraint  should  be  in  the  initial  working  set  at  its  lower  bound. 

2  The  constraint  should  be  in  the  initial  working  set  at  its  upper  bound. 

3  The  constraint  should  be  in  the  initial  working  set  as  an  equality.  This 

value  must  not  be  specified  unless  BL(j')  =  BU(j).  The  values  I,  2  or  3 
all  have  the  same  cITcct  when  BL(j)  =  BU(j). 
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Other  values  of  ISTATE  are  also  acceptable.  In  particular,  if  QPSOL  has  been  called 
previously  with  the  same  values  of  N  and  NCLIN,  ISTATE  already  contains  satisfactory 
values. 

When  QPSOL  exits  with  INFORM  set  to  0,  1  or  3,  the  values  in  the  array  ISTATE  indicate 
the  status  of  the  constraints  in  the  active  set  at  the  solution.  Otherwise,  ISTATE 
indicates  the  composition  of  the  working  set  at  the  final  iterate. 

X  is  a  real  array  of  dimension  N  that  contains  the  current  estimate  of  the  solution.  On 

entry  to  QPSOL,  X  must  be  defined;  on  exit  from  QPSOL,  X  contains  the  best  estimate  of 
the  solution. 
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INFORM  is  an  integer  that  indicates  the  result  of  QPSOL.  (When  MSGLVL  >  0,  a  short  description 
of  INFORM  is  printed.)  The  possible  values  of  INFORM  are: 

INFORM  Definition 

0  X  is  a  strong  local  minimum,  i.c.,  the  projected  gradient  is  negligible, 
the  Lagrange  multipliers  are  optimal,  and  the  projected  llcssian  is 
positive  definite.  In  some  cases,  a  zero  value  of  INFORM  means  that  X  is 
a  global  minimum  (c.g.,  when  the  llcssian  matrix  is  positive  definite). 

1  X  is  a  weak  local  minimum  (the  projected  gradient  is  negligible,  the 
Lagrange  multipliers  are  optimal,  but  the  projected  llcssian  is  only 
positive  semi-definite).  This  means  that  the  solution  is  not  unique. 

2  The  solution  appears  to  be  unbounded,  i.c.,  the  quadratic  function  is 
unbounded  below  in  the  feasible  region.  This  value  of  INFORM  occurs 
when  a  step  of  infinity  would  have  to  be  taken  in  order  to  continue  the 
algorithm. 

3  X  appears  to  be  a  local  minimum,  but  optimality  cannot  be  verified 
because  some  of  the  Lagrange  multi  pliers  are  very  small  in  magnitude. 

4  The  iterates  of  the  QP  phase  could  be  cycling,  since  a  total  of  50 
changes  were  made  to  the  working  set  without  altering  X. 

5  The  limit  of  ITMAX  iterations  was  reached  in  the  QP  phase  before 
normal  termination  occurred. 

6  The  LP  phase  terminated  without  finding  a  feasible  point,  ami  hence 
it  is  not  possible  to  satisfy  all  the  constraints  to  within  the  tolerances 
specified  by  the  FEATOL  array.  In  this  ease,  the  final  iterate  will  reveal 
values  for  which  there  will  be  a  feasible  point  (c.g.,  a  feasible  point  will 
exist  if  the  feasibility  tolerance  for  each  violated  constraint  exceeds 
its  RESIDUAL  at  the  final  point).  The  modified  problem  (with  altered 
values  in  FEATOL)  may  then  be  solved  using  a  warm  start. 

7  The  iterates  may  be  cycling  during  the  LP  phase;  see  the  comments 
above  under  INFORM  =  4. 

8  The  limit  of  ITMAX  iterations  was  reached  during  the  LP  phase. 

9  An  input  parameter  is  invalid. 

ITER  is  an  integer  that  gives  the  number  of  iterations  performed  in  either  the  LP  phase  or 
the  QP  phase,  whichever  was  last  entered.  (Note  that  ITER  is  reset  to  zero  after  the 
LP  phase.) 
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OBJ 

CLAMDA 


is  the  value  of  the  quadratic  objective  function  at  X  if  X  is  feasible  (INFORM  <  5),  or  the 
sum  of  infeasibilities  at  X  otherwise  (6  <  INFORM  <  8). 

is  a  real  array  of  dimension  NCTOTL  that  contains  the  Lagrange  multiplier  for  every 
constraint  with  respect  to  the  current  working  set.  The  ordering  of  CLAMDA  follows  the 
convention  given  above  under  BL,  i.e.,  the  first  N  components  contain  the  multipliers 
for  the  bound  constraints  on  the  variables,  and  the  remaining  components  contain  the 
multipliers  for  the  general  linear  constraints.  If  ISTATE(jf)  =  0  (i.e.,  constraint  j  is  not 
in  the  working  set),  CLAMDA(j)  is  zero.  If  X  is  optimal,  CLAMDA^)  should  be  non-negative 
if  ISTATE(j)  =  1  and  non-positive  if  ISTATE(jf)  =  2. 
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IV  is  an  integer  array  of  dimension  LENIV,  which  provides  integer  workspace  Tor  QPSOL. 
LENIW  is  the  dimension  of  IV,  and  must  be  at  least  N  +  2  +  min  (N,  NCLIN). 

V  is  a  real  array  of  dimension  LENV,  which  provides  real  workspace  for  QPSOL. 

LENV  is  the  dimension  of  V.  If  LP  ==  .FALSE,  or  NCLIN  >  N,  LENV  must  be  at  least  2N*  +  4N  + 

NROVA+  2NC0N,  where  NCON  =  max(t, NCLIN).  If  LP  =  .TRUE,  and  NCLIN  <  N,  LENV 
must  be  at  least  2NC0N2  +  -IN  +  NROVA  +  2 NCON. 

IfMSGLVL  >  0,  the  amounts  of  workspace  provided  and  required  arc  printed.  As  an  alternative 
to  eompuling  LENV  from  the  formula  given  above,  the  user  may  prefer  to  obtain  an  appropriate 
value  from  the  output  of  a  preliminary  run  with  a  positive  value  of  MSGLVL  and  LENV  set  to  1 
(QPSOL  will  then  terminate  with  INFORM  =  9). 
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The  subroutines  associated  specifically  with  the  QPSOL  package  arc  the  following: 


ADDCON 

ALLOC 

BDPERT 

BNDALF 

CHKDAT 

DELCON 

FINDP 

GETLAM 

LPBGST 

LPCORE 

LPCRSH 

LPDUMP 

LPGRAD 

LPPRT 

MOVEX 

QPCHKP 

QPCOLR 

QPCORE 

QPCRSH 

QPDUMP 

QPGRAD 

QPPRT 

PRTSOL 

RSOLVE 

TQADD 

TSOLVE 

ZYPROD. 

QPSOL  also  uses  the  basic  linear  algebra  subroutines 


AXPY 

CONDVC 

COPYMX 

COPYVC 

DOT 

DSCALE 

ELM 

ELMGEN 

ETAGEN 

QUOTNT 

REFGEN 

R0T3 

ROTGEN 

SSCALE 

V2N0RM 

ZEROVC 

and  the  subroutine  MCHPAR,  which  defines  machine-dependent  constants  (sec  Section  11). 

The  subroutines  in  the  QPSOL  package  use  the  following  labelled  COMMON  areas: 

SOLMCH  (15  REAL  variables;  sec  Section  11) 

S0L1CM  (3  INTEGER  variables) 

S0L3CM  (1  INTEGER  variables) 

S0L4CM  (10  REAL  variables) 

S0L5CM  (3  REAL  variables) 

S0L1LP  (15  INTEGER  variables) 

S0L2LP  (l  LOGICAL  variable.) 
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•.  DESCRIPTION  OF  THE  PRINTED  OUTPUT 


This  section  describes  the  intermediate  printout  produced  by  QPSOL.  When  MSGLVL  >  5,  a 
line  or  output  is  produced  for  every  change  in  the  working  set  (thus,  several  lines  may  bo  printed 
during  a  single  iteration). 

To  aid  interpretation  of  the  printed  results,  we  mention  the  convention  for  numbering  the 
constraints:  indices  1  through  N  refer  to  the  bounds  on  the  variables,  and  indices  N  +  1  through 
N  4-  NCLIN  refer  to  the  general  constraints.  When  the  status  of  a  constraint  changes,  the  index 
of  the  constraint  is  printed,  along  with  the  designation  “L”  (lower  bound),  “U”  (upper  bound)  or 
“E”  (equality),  ir  the  problem  is  non-convcx,  the  character  “V”  may  appear  alongside  an  index 
in  the  “delete”  column.  This  will  occur  if  the  initial  projected  Hessian  is  not  sufficiently  positive 
definite  (and  therefore  the  Cholcsky  factor  corresponds  only  to  a  subset  of  the  columns  of  Z;  see 
Section  2).  The  “V”  is  used  to  indicate  that  the  Cholcsky  factor  has  been  expanded  to  include  a 
new  column  or  Z.  The  associated  index  gives  the  current  dimension  or  the  Cholcsky  factor. 

In  the  LP  phase,  the  printout  includes  the  following: 


ITU 

is  the  iteration  count. 

KDEL 

is  the  index  of  the  constraint  deleted  from  the  working  set.  If  KDEL  is 
zero,  no  constraint  was  deleted. 

KADD 

is  the  index  or  the  constraint  added  to  the  working  set.  If  KADD  is  zero, 
no  constraint  was  added. 

STEP 

is  the  step  taken  along  the  computed  search  direction. 

MUMIMF 

is  the  number  or  violated  constraints  (infcasibilitics). 

SUM IMF 

is  a  weighted  sum  of  the  magnitudes  or  the  constraint  violations. 

LPOBJ 

is  the  value  of  the  linear  objective  function  erx.  It  is  printed  only  if  LP 
is  .TRUE. 

• 

During  the  QP  phase,  the  printout  includes  the  following: 

ITU 

is  the  iteration  count  (reset  to  zero  after  the  LP  phase). 

KDEL 

is  the  index  of  the  constraint  deleted  from  the  working  set.  If  KDEL  is 
zero,  no  constraint  was  deleted. 

KADD 


is  the  index  of  the  constraint  added  to  the  working  set.  If  KADD  is  zero, 
no  constraint  was  added. 
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STEP 


NHES8 

OBJECTIVE 

NCOLZ 


is  the  step  aic  taken  along  the  direction  of  search  (if  STEP  is  1.0,  the 
current  point  is  a  minimum  in  the  subspacc  defined  by  the  current 
working  set). 

is  the  number  of  calls  to  subroutine  QPHESS. 
is  the  value  of  the  quadratic  objective  function. 

is  the  number  of  columns  of  Z  (see  Section  2).  In  general,  it  is  the 
dimension  of  the  subspace  in  which  the  quadratic  is  currently  being 
minimized. 


NORM  GFREE  is  the  Euclidean  norm  of  the  gradient  of  the  objective  function  with 

respect  to  the  free  variables,  i.c.  variables  not  currently  held  at  a  bound 
(NORM  GFREE  is  not  printed  if  ORTHOG  is  .FALSE.).  In  some  cases,  the 
objective  function  and  gradient  arc  updated  raLhcr  than  recomputed.  If 
so,  this  entry  will  be  “ — ”  to  indicate  that  the  gradient  with  respect  to 
the  free  variables  has  not  been  computed. 


NORM  QTG 


NORM  ZTG 
HESS  MOD 


is  a  weighted  norm  of  the  gradient  of  the  objective  function  with  respect 
to  the  free  variables  (NORM  QTG  is  not  printed  ir  ORTHOG  is  .TRUE.).  In 
some  cases,  the  objective  function  and  gradient  arc  updated  rather  than 
recomputed.  If  so,  this  entry  will  be  “ — ”  to  indicate  that  the  gradient 
with  respect  to  the  free  variables  has  not  been  computed. 

is  the  Euclidean  norm  or  the  projected  gradient  (sec  Section  2). 

is  the  correction  added  to  the  diagonal  of  the  projected  Hessian  to  ensure 
that  a  satisfactory  Cholcsky  factorization  exists  (sec  Section  2).  When 
the  projected  Hessian  is  sufficiently  positive  definite,  HESS  MOD  will  be 
zero. 


When  MSGLVL  =  1  or  MSGLVL  >  10,  the  summary  printout  at  the  end  of  execution  of  QPSOL 
includes  a  listing  of  the  status  of  every  constraint.  Note  that  default  names  arc  assigned  to  all 
variables  and  constraints. 

The  following  describes  the  printout  for  each  variable. 


VARIABLE  is  the  name  (VARBL)  and  index  j  or  the  variable. 

STATE  gives  the  state  of  the  variable  (FR  ir  neither  bound  is  in  the  working  set, 

EQ  if  a  fixed  variable,  LL  if  on  its  lower  bound,  UL  if  on  its  upper  bound). 
If  VALUE  lies  outside  the  upper  or  lower  bounds  by  more  than  FEATOL(j), 
STATE  will  be  “++”  or  respectively. 


! 
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VALUE 

LOVER  BOUND 

UPPER  BOUND 
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is  the  value  of  the  variable  at  the  final  iteration. 

is  the  lower  bound  specified  for  the  variable.  (“NONE”  indicates  that 
BL(j)  <  — BIGBND.) 

is  the  upper  bound  specified  for  the  variable.  (“NONE”  indicates  that 
BU(j')  >  BIGBND.) 


LAGR  MULTIPLIER  is  the  value  of  the  Lagrange  multiplier  Tor  the  associated  bound  con¬ 
straint.  This  will  be  zero  if  STATE  is  FR.  If  X  is  optimal,  the  multiplier 
should  be  non-negative  if  STATE  is  LL,  and  non-positive  if  STATE  is  UL. 

RESIDUAL  is  the  difference  between  the  variable  and  the  nearer  of  its  bounds  BL (j) 

and  BU(j). 


The  following  summary  printout  is  given  for  each  general  constraint. 


LINEAR  CONSTR 
STATE 


VALUE 

LOWER  BOUND 

UPPER  BOUND 

LAGR  MULTIPLIER 

RESIDUAL 


is  the  name  (LNCON)  and  index  i,  t  =  1  to  NCLIN,  of  the  constraint. 

is  the  state  of  the  constraint  (FR  for  a  constraint  not  in  the  working  set, 
EQ  for  an  equality,  LL  for  an  inequality  constraint  at  its  lower  bound,  UL 
for  an  inequality  constraint  at  its  upper  bound).  If  VALUE  lies  outside 
the  upper  or  lower  bounds  by  more  than  its  feasibility  tolerance,  STATE 
will  be  or  " — ”  respectively. 

is  the  value  of  the  constraint  at  the  final  point,  i.c.,  the  appropriate 
component  of  the  vector  Ax. 

is  the  specified  lower  bound  for  the  constraint.  (“NON E"  indicates  that 
BL(N  +  t)  <  -BIGBND.) 

is  the  specified  upper  bound  for  the  constraint.  (“NONE”  indicates  that 
BU(N  +  t)  >  BIGBND.) 

is  the  value  of  the  Lagrange  multiplier.  This  will  bo  zero  if  STATE  is  FR. 
IfX  is  optimal,  the  multiplier  should  be  non-negative  ir  STATE  is  LL,  and 
non-positive  if  STATE  is  UL. 

is  the  residual  of  the  constraint  with  respect  to  its  nearer  bound,  i.e., 
the  difference  between  VALUE  and  the  nearer  of  its  two  bounds. 


3  ^ 


>5 
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Reason  for  termination 


Recommended  Action 


Underflow 


Overflow 


INFORM  =  3 


INFORM  =  A 


If  the  machine  parameter  indicating  an  underflow  check  (WMACH(9))  is 
zero,  floating-point  underflow  may  occur  occasionally,  but  can  usually  be 
ignored.  To  avoid  underflow,  set  WMACH(O)  to  a  positive  value;  however, 
this  will  lead  to  a  noticeable  loss  of  eflicicncy.  If  underflow  continues  to 
occur  for  no  apparent  reason,  contact  the  authors  at  Stanford  University. 

If  the  printed  output  before  the  overflow  error  contains  a  warning  about 
serious  ill-conditioning  in  the  working  set  when  adding  the  j'-th  con¬ 
straint,  it  may  be  possible  to  avoid  the  difficulty  by  increasing  the  mag¬ 
nitude  of  FEATOL (j)  and  rerunning  the  program.  If  the  message  recurs 
even  after  this  change,  the  olTcnding  linearly  dependent  constraint  (with 
index  uj”)  must  be  removed  from  the  problem.  If  a  warning  message 
did  not  precede  the  fatal  overflow,  contact  the  authors  at  Stanford 
University. 

QPSOL  has  probably  found  a  solution.  However,  the  presence  of  very 
small  Lagrange  multipliers  means  that  the  predicted  active  set  may  be 
incorrect,  or  that  X  may  be  only  a  constrained  stationary  point  rather 
than  a  local  minimum.  The  method  in  QPSOL  is  not  guaranteed  to 
find  the  correct  active  set  when  there  arc  small  multipliers.  QPSOL 
attempts  to  delete  constraints  with  zero  multipliers,  but  this  docs  not 
necessarily  resolve  the  issue.  The  determination  of  the  correct  active  set 
is  a  combinatorial  problem  that  may  require  an  extremely  large  amount 
of  lime.  The  occurrence  of  small  multipliers  often  ( but  not  always ) 
indicates  that  there  arc  redundant  constraints. 

This  value  will  occur  if  50  iterations  arc  performed  in  the  QP  phase 
without  changing  X.  The  user  should  check  the  printed  output  for  a 
repeated  pattern  of  constraint  deletions  and  additions,  ir  a  sequence  of 
constraint  changes  is  being  repeated,  the  iterates  arc  probably  cycling. 
(QPSOL  docs  not  contain  a  method  that  is  guaranteed  to  avoid  cycling, 
which  would  be  combinatorial  in  nature.)  Cycling  may  occur  in  two 
circumstances:  at  a  constrained  stationary  point  where  there  arc  some 
small  or  zero  Lagrange  multipliers  (sec  the  discussion  of  INFORM  =  3); 
or  at  a  point  (usually  a  vertex)  where  the  constraints  that  arc  satisfied 
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INFORM  =  5 


INFORM  =  6 


exactly  are  nearly  linearly  dependent.  In  the  latter  case,  the  user  has  the 
option  of  identifying  the  offending  dependent  constraints  and  removing 
them  from  the  problem,  or  restarting  the  run  with  larger  values  of 
FEATOL  for  nearly  dependent  constraints.  If  QPSOL  terminates  with 
INFORM  =  4,  but  no  suspicious  pattern  of  constraint  changes  can  be 
observed,  it  may  be  worthwhile  to  restart  with  the  final  X  (with  or 
without  the  warm  start  option). 

The  value  of  ITMAX  may  be  too  small.  If  the  method  appears  to  be  mak¬ 
ing  progress  (c.g.,  the  objective  function  is  being  satisfactorily  reduced), 
increase  ITMAX  and  rerun  QPSOL  (possibly  using  the  warm  start  facility 
to  specify  the  initial  working  set).  If  ITMAX  is  already  large,  but  some  of 
the  constraints  could  be  nearly  linearly  dependent,  check  the  output  for 
a  repeated  pattern  of  constraints  entering  and  leaving  the  working  set. 
(Near-dependencies  arc  often  indicated  by  wide  variations  in  size  in  the 
diagonal  elements  of  the  T  matrix,  which  will  be  printed  if  MSGLVL  > 
30.)  In  this  case,  the  algorithm  could  be  cycling  (sec  the  comments  for 
INFORM  =  4.) 

The  LI’  phase  has  terminated  without  finding  a  feasible  point,  which 
means  that  no  feasible  point  exists  for  the  given  FEATOL  array.  The  user 
should  check  that  there  arc  no  constraint  redundancies.  If  the  data  for 
the  j-th  constraint  arc  accurate  only  to  the  absolute  precision  S,  the 
user  should  ensure  that  the  value  of  FEATOL(j)  is  greater  than  6.  For 
example,  if  all  elements  or  A  arc  or  order  unity  and  arc  accurate  only 
to  three  decimal  places,  every  component  of  FEATOL  should  be  at  least 
10~3. 


INF0RM=  7  or  8 


These  values  arc  the  analogue  in  the  LP  phase  procedure  or  INFORM 
values  4  and  5. 
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This  program  has  been  written  in  ANSI  (1966)  Fortran  and  tested  on  an  IBM  3081  computer 
using  the  WATFIV  Compiler  Version  1  Level  6.  All  subroutines  in  QPSOL  are  PFORT-compatible 
(Ryder,  1974),  except  Tor  CHKDAT,  GETLAM  and  PRTSOL,  which  contain  A2  Format  specifications. 

At  the  beginning  of  QPSOL,  the  subprogram  MCHPAR  is  called  to  assign  various  machine- 
dependent  parameters.  These  parameters  arc  stored  in  the  array  WMACH(15)  in  the  labelled  COMMON 
block  SOLMCH. 

The  specification  of  MCHPAR  is 


SUBROUTINE  MCHPAR 

REAL  WMACH 

COMMON  /SOLMCH/  WMACH (IS) 


The  first  eleven  components  of  the  REAL  array  WMACH  must  be  set  in  MCHPAR.  The  components 
of  WMACH  arc  defined  as  follows. 


WMACH(I) 

WMACH(2) 

WMACH(3) 

WMACH(4) 

WMACH(5) 

WMACH(6) 

WMACH(7) 

WMACH(8) 

WMACH(9) 


WMACH(IO) 

WMACH(ll) 


Definition 

is  NBASE,  the  base  of  floating-point  arithmetic. 

is  NDIGIT,  the  number  of  NBASE  digits  of  precision. 

is  EPSMCH,  the  floating-point  precision. 

is  RTEPS,  the  square  root  or  EPSMCH. 

is  FLMIN,  the  smallest  positive  Moating- point  number. 

is  RTMIN,  the  square  root  of  FLMIN. 

is  FLMAX,  the  largest  positive  Moating-point  number. 

is  RTMAX,  the  square  root  of  FLMAX. 

is  UNDFLW,  which  specifies  whether  or  not  NPSOL  should  check  for 
underflow  in  certain  computations.  If  UNDFLW  =  0,  no  underflow 
checking  will  be  performed.  If  UNDFLW  is  set  to  a  positive  number, 
QPSOL  will  check  for  underflow  and  will  replace  too-small  quantities 
by  r.cro.  Note  that  QPSOL  will  run  faster  if  no  underflow  checking 
takes  place ,  i.c.  if  W1IACH(9)  =  0.0. 

is  NIN,  the  file  number  Tor  the  input  stream. 

is  NOUT,  the  file  number  for  the  output  stream. 
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The  following  version  of MCHPAR  (which  is  provided  hy  the  Systems  Optimization  Laboratory) 
contains  the  parameters  associated  with  double  precision  on  a  machine  in  the  II5M  i570  series. 
Tin)  user  must  substitute1  a  version  of  MCHPAR  Lh.it  is  appropriate  /‘or  the  machine  to  he  used. 


SUBROUTINE  MCHPAR 
C 

DOUBLE  PRECISION  WMACH 
COMMON  /SOLMCH/  WMACH ( IS) 

C 

C  MCHPAR  MUST  DEFINE  THE  RELEVANT  MACHINE  PARAMETERS  AS  FOLLOWS. 

C  WMACH! 1  )  =  NBASE  =  BASE  OF  FLOATING-POINT  ARITHMETIC. 

C  WMACH! 2 )  =  NDIGIT  =  NO.  OF  BASE  WMACH! 11  DIGITS  OF  PRECISION. 

C  WMACH! 3)  =  EPSMCH  =  FLOATING-POINT  PRECISION. 

C  WMACH! 4)  =  RTEPS  =  SORT!  EPSMCH I . 

C  WMACH! 5)  =  FLMIN  =  SMALLEST  POSITIVE  FLOATING-POINT  NUMBER. 

C  WMACH! 6)  =  RTMIN  =  SORT! FLMIN). 

C  WMACH! 7)  =  FLMAX  =  LARGEST  POSITIVE  FLOATING-POINT  NUMBER. 

C  WMACH! 8)  =  RTMAX  =  SORT! FLMAX). 

C  WMACH! 9)  =  UNDFLW  =  0.0  IF  UNDERFLOW  IS  NOT  FATAL,  *VE  OTHERWISE. 

C  WMACH! 10)  =  NIN  =  STANDARD  FILE  NUMBER  OF  THE  INPUT  STREAM. 

C  WMACH! 11)  =  NOUT  =  STANDARD  FILE  NUMBER  OF  THE  OUTPUT  STREAM. 

C 

INTEGER  NBASE,  NDIGIT.  NIN,  NOUT 

DOUBLE  PRECISION  OSQRT 
C 

NBASE  =  16 

NDIGIT  =  14 

WMACH! 1 )  =  NBASE 

WMACH! 2)  =  NDIGIT 

WMACH! 3  I  =  WMACH! 11**11  -  NDIGIT) 

WMACH! 4 )  =  OSQRT! WMACH! 3 ) I 

WMACH! 5 )  =  WMACH! I )»*! -62) 

WMACH! 6)  =  DSQRT! WMACH! 5  )  ) 

WMACH 1 7 )  =  WMACH! 1 )**61 

WMACH! 8)  =  OSQRT! WMACH! 7)) 

WMACH! 9)  =  0.00*0 

NIN  =  5 

NOUT  =  6 

WMACH! 10)  =  NIN 
WMACH!  11)=  NOUT 
C 

C -  IN  WATFIV,  ALLOW  UP  TO  100  UNDERFLOWS. 

C -  CALL  TRAPS  !  0,0,100  ) 

RETURN 

C 

C  END  OF  MCHPAR 
END 
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The  values  of  NBASE,  HDIGIT,  EPSMCH,  FLMIN  and  FLMAX  for  several  machines  arc  given  in  the 
following  tabic,  for  both  single  and  double  precision;  RTEPS,  RTUIN  and  RTMAX  may  be  computed 
using  Fortran  statements.  The  values  NIN  and  NOUT  depend  on  the  machine  installation. 

For  each  precision,  we  give  two  values  for  EPSMCH,  FLMIN  and  FLMAX.  The  first  value  is  a 
Fortran  decimal  approximation  of  the  exact  quantity;  use  or  this  value  in  MCHPAR  should  cause 
no  difficulty  except  in  extreme  circumstances.  The  second  value  is  the  exact  mathematical 
representation. 


NBASE 


NDIGIT 


NBASE 


NDIGIT 


Table  of  machine-dependent  parameters 


Variable  IBM  360/370  CDC  6000/7000 


Single 


Single 


DEC  10/20 
Single 


Uni  vac  1100 
Single 


DEC  VAX 
Single 


2 


Variable  IBM  360/370  CDC  6000/7000 


Double 


Double 


DEC  10/20 
Double 


Univac  1100 
Double 


DEC  VAX 
Double 


*1  _ 
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12.  EXAMPLE  PROGRAM  AND  OUTPUT 

This  section  contains  a  listing  and  the  computed  results  from  a  sample  main  program  that 
calls  QPSOL  to  solve  an  indefinite  quadratic  program.  The  problem  has  seven  variables  and  seven 
general  constraints. 

The  vector  c  is  given  by 


c  =  (-.02,  -.2,  -.2,  -.2,  -.2,  .04,  .04)r. 


The  Hessian  is 


0  0  0  0  0  -2  -2 

V  0  0  0  0  0  -2  -2  / 

and  is  defined  by  the  subroutine  QPHES1,  which  docs  not  store  //  explicitly. 
The  general  constraint  matrix  A  is 


1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

.15 

.04 

.02 

.04 

.02 

.01 

.03 

.03 

.05 

.08 

.02 

.06 

.01 

0.0 

.02 

.04 

.01 

.02 

.02 

0.0 

0.0 

.02 

.03 

0.0 

0.0 

.01 

0.0 

0.0 

.70 

.75 

.80 

.75 

.80 

.97 

0.0 

.02 

.06 

.08 

.12 

.02 

.01 

.97 

tors  l  and  u  are 

•l,  - 

-.01, 

-.04, 

-•1, 

-.01, 

-.01, 

-oo, 

-00, 

-00, 

-oo, 

1 

of 

s 

1 

.003] 

.03, 

.02, 

.05,  +oo,  +oo, 

31  — ■ 


-.13,  -.0049,  -.0064,  -.0037,  -.0012,  +oo,  .002  )r. 
The  starting  point  xq  (which  is  infeasible)  is 

x0  =  (-.01,  -.03,  0.0,  -.01,  -.1,  .02,  .01)r. 


,T  «  *,* 
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The  computed  solution  (to  live  figures)  is 

x  —  (-.01,  -.069865,  .018259,  -.024281,  -.062006,  .013805,  .0040665  )T. 
One  bound  constraint  and  four  general  constraints  are  active  at  the  solution. 


s  v 
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c»»*k*begin  file  qpmain  fortran  d. 
c 

C  EXAMPLE  PROGRAM  FOR  SUBROUTINE  QPSOL. 

C  DOUBLE  PRECISION  VERSION  3.2.  SEPTEMBER  1984. 

C  THE  VALUE  OF  THE  PARAMETER  FEATOL  IS  APPROPRIATE  FOR  A  MACHINE 


C 

WITH  A  PRECISION  OF 

15  DECIMAL  DIGITS. 

C 

******* *  *  *****************************  **  ************************ ***** 

mrn 

U  . 

1 

INTEGER 

I,  INFORM,  ITER,  ITMAX,  J,  LIWORK,  LWORK 

2 

INTEGER 

MSGLVL,  N,  NCLIN,  NCOLH,  NCOLH1 ,  NCTOTL 

•  * 

3 

INTEGER 

NIN,  NOUT,  NROWA,  NROWH,  NROWH 1 

4 

INTEGER 

ISTATE( 14),  IWORKt  50 ) 

5 

DOUBLE  PRECISION 

BIGBND,  EPSMCH,  OBJ,  RTEPS 

6 

DOUBLE  PRECISION 

ZERO,  TWO 

7 

DOUBLE  PRECISION 

A( 7,7) ,  BL(  14),  BU(  14) ,  CLAMDA(  14 ) ,  CVEC(7) 

8 

OOUBLE  PRECISION 

FEATOLf 14 ),  HESS( 1,1),  HESS1 ( 7,7) ,  X(7> 

9 

DOUBLE  PRECISION 

WORM  200) 

10 

DOUBLE  PRECISION 

DSQRT  1 

11 

LOGICAL 

COLD,  LP,  ORTHOG 

12 

EXTERNAL 

QPHES1 ,  QPHES2 

•  .* 

13 

DATA 

ZERO  ,  TWO 

* 

/0.00+0,  2.0D+0/ 

c 

►i'- 

c 

SET  THE  OECLARED  ARRAY  DIMENSIONS. 

m 

c 

NROWA  =  THE  DECLARED  ROW  DIMENSION  OF  A. 

c 

NROWH  *  THE  DECLARED  ROW  DIMENSION  OF  HESS. 

.. 

c 

NCOLH  -  THE  NUMBER 

OF  COLUMNS  IN  HESS. 

*,* 

c 

(IF  QPHESS  OEALS  WITH  THE  HESSIAN  IMPLICITLY, 

c 

NROWH  AND 

NCOLH  CAN  BOTH  BE  1 . ) 

c 

LIWORK  =  THE  LENGTH 

OF  THE  INTEGER  WORK  ARRAY. 

•T» 

c 

LWORK  -  THE  LENGTH 

OF  THE  OOUBLE  PRECISION  WORK  ARRAY. 

c 

S 

r 

14 

NROWA  -  7 

15 

NROWH  =  1 

fJ- 

16 

NCOLH  =  1 

17 

LIWORK  =  50 

18 

r 

LWORK  =  200 

V* 

.-.T, 

c 

SET  THE  APPROXIMATE 

MACHINE  PRECISION. 

c 

■ 

19 

r 

EPSMCH  =  1.0D-15 

y! 

.*** 

L 

c 

ALLOW  UP  TO  20  ITERATIONS  TO  FIND  A  FEASIBLE  POINT, 

C 

r 

AND  THE  SAME  NUMBER 

TO  MINIMIZE  THE  QUADRATIC  FUNCTION. 

•/ 

20 

ITMAX  =  20 

1  > 

C 

W  ■ 

c 

ASK  FOR  BRIEF  OUTPUT  EACH  ITERATION,  AND  A  FULL  PRINT-OUT 

c 

r 

OF  THE  FINAL  SOLUTION. 

21 

L 

MSGLVL  =  10 

•V 

SET  THE  PROBLEM  DIMENSIONS. 

N  =  THE  NUMBER  OF  VARIABLES. 

NCLIN  =  THE  NUMBER  OF  GENERAL  LINEAR  CONSTRAINTS  (MAY  BE  0). 
NCTQTL  =  THE  TOTAL  NUMBER  OF  VARIABLES  AND  GENERAL  CONSTRAINTS. 

(THE  ARRAYS  ISTATE.  BL,  BU,  C LAMBDA  MUST  BE  AT  LEAST 
THIS  LONG.) 


V  v, 

V, 
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N  =  7 

NCLIN  =  7 
NCTOTL  =  N  ♦  NCLIN 

BOUNDS  GREATER  THAN  BI68ND  HILL  BE  TREATED  AS  PLUS  INFINITY. 

BOUNDS  LESS  THAN  -  BI6BND  MILL  BE  TREATED  AS  MINUS  INFINITY. 

BI6BND  =  1 . 0E*10 

ANY  BOUND  OR  LINEAR  CONSTRAINT  NAY  BE  VIOLATED  BY  AS  MUCH  AS  FEATOL. 

RTEPS  =  DSQRT(  EPSNCH  ) 

DO  20  J  =  1.  NCTOTL 
FEATOL(J)  =  RTEPS 
20  CONTINUE 

A  COLO  START  IS  NEEDED  FOR  THE  FIRST  CALL  TO  QPSOL. 

HE  MANY  TO  SOLVE  A  QUADRATIC  PROGRAM.  NOT  AN  LP  PROBLEM. 

USE  AN  ORTHOGONAL  FACTORIZATION  OF  THE  MATRIX  OF  CONSTRAINTS 
IN  THE  N0RKIN6  SET. 


30 

COLD  =  .TRUE. 

31 

LP  =  .FALSE 

32 

ORTHOG  *  .TRUE. 

READ  THE  DATA  ARRAYS. 


THE  UNIT  NUTSER  FOR  INPUT. 

THE  UNIT  NUMBER  FOR  PRINTING. 

THE  LINEAR  PART  OF  THE  OBJECTIVE  FUNCTION. 
THE  6ENERAL  CONSTRAINT  MATRIX. 

THE  LOWER  BOUNDS  ON  X  AND  A«X. 

THE  UPPER  BOUNDS  ON  X  AND  A*X. 

THE  INITIAL  ESTIMATE  OF  THE  SOLUTION. 


33 

NIN 

=  5 

34 

NOUT 

=  6 

35 

READ 

(NIN. 

1000) 

( 

CVEC(J), 

J=1  ,N  ) 

36 

READ 

(NIN, 

1000) 

(( 

A(I,J), 

J=1,N  ),  1=1 

37 

READ 

(NIN, 

1000) 

( 

BL1J), 

J=1, NCTOTL  ) 

36 

READ 

(NIN, 

1000) 

( 

BU( J ) , 

J=1, NCTOTL) 

39 

READ 

(NIN, 

1000) 

( 

X(J), 

J=1  ,N) 

•. ♦ 


PRINT  THE  DATA. 

IF  (NOUT  .LE.  0)  60  TO  50 

WRITE  (NOUT,  2000)  (CVEC(J),  J=1,N) 

WRITE  (NOUT,  2100)  ((A(I,J),  J=1,N),  1=1, NCLIN) 
WRITE  (NOUT,  2200)  (  BL(J),  J=1 , NCTOTL) 

WRITE  (NOUT,  2300)  (  BU(J),  J=1 , NCTOTL) 

WRITE  (NOUT,  2400)  (  X(J),  J=1,N) 


SOLVE  THE  PROBLEM. 

THE  HESSIAN  IS  DEFINED  IMPLICITLY  BY  SUBROUTINE  QPHES1 . 

50  CALL  QPSOL1  ITMAX,  MS6LVL,  N, 

•  NCLIN.  NCTOTL,  NROMA,  NROWH,  NCOLH, 

•  BI6BND,  A,  BL,  BU,  CVEC,  FEATOL,  HESS,  RPHES1 , 

•  COLD,  LP,  ORTHOG,  1ST ATE,  X, 

•  INFORM,  ITER,  OBJ,  CLAMDA, 

•  I WORK,  LI  WORK,  WORK,  LMORK  ) 


tz 

-  .1 


V.'-o 

V.- 


'j 


Ql‘  'SOI./38 


12.  KXAMPI.K  I'KOGIIAM  AND  OUTPUT 


C 

C  TEST  FOR  AN  ERROR  CONDITION. 

C 

47  IF  (INFORM  .6T.  0)  60  TO  900 

C 
C 

C  THE  FOLLOWING  IS  FOR  ILLUSTRATIVE  PURPOSES  ONLY. 

C  HE  DO  A  HARM  START  NITH  THE  FINAL  WORKING  SET  OF  THE  PREVIOUS  RUN. 
C  THIS  TIME  HE  STORE  THE  HESSIAN  EXPLICITLY  IN  HESSt  > 

C  AND  USE  THE  CORRESPONDING  SUBROUTINE  QPHES2. 

C 


40 

HRITE  (NOUT,  2500) 

49 

COLD  ■  .FALSE. 

50 

MSGLVL  =  5 

51 

NROHH t  s  7 

52 

C 

NCOLHt  =  7 

53 

DO  200  J  :  1,  N 

54 

DO  100  I  i  1,  N 

55 

HESSKI.J)  *  ZERO 

56 

100 

CONTINUE 

57 

IF  (J  ,LE.  5)  HESSKJ.J) 

50 

IF  (J  .6T.  5)  HESSKJ.J) 

59 

C 

200 

CONTINUE 

60 

HESS1(3,4)  s  TWO 

61 

HESSK4.3)  *  TWO 

62 

HESS1 (6,7)  *  -  TWO 

63 

C 

HESS1(7,6)  *  -  TWO 

64 

CALL  QPSOL!  ITMAX,  MSGLVL,  1 

•  NCLIN»  NCTOTL,  NR OH A,  NR0NH1 ,  NCOLHt , 

•  BIGBND,  A»  BL.  BU,  CVEC,  FEATOL,  HESS1 ,  QPHES2, 

•  COLD.  LP.  ORTHOG.  ISTATE.  X. 

•  INFORM,  ITER,  OBJ,  CLAMDA, 

•  I WORK,  LINORK,  WORK,  LHORK  ) 

C 

65  IF  (INFORM  .6T.  0)  GO  TO  900 

66  STOP 
C 

C  ERROR  EXIT. 

C 

67  900  WRITE  (NOUT,  3000)  INFORM 

60  STOP 

C 

69  1000  FORMAT( 7E10.2 ) 

70  2000  FORMAT( /  14H  CVEC.  /  (IX,  7F1Q.2>> 

71  2100  FORMAT! /  14H  ROWS  OF  A.  /  (IX,  7F10.2M 

72  2200  FORMAT) /  14H  LOWER  BOUNDS.  /  (IX,  7EI0.2)) 

73  2300  FORMAT!/  14H  UPPER  BOUNDS.  /  (IX,  7E10.2)) 

74  2400  FORMAT!/  12H  INITIAL  X.  /  (IX,  7F10.2)) 

75  2500  FORMAT! //46H  A  RUN  OF  THE  SAME  EXAMPLE  WITH  A  HARM  START....) 

76  3000  FORMAT) /  32H  QPSOL  TERMINATED  WITH  INFORM  *,  13) 

C 

C  END  OF  THE  EXAMPLE  PROGRAM  FOR  QPSOL. 

77  END 


70 

79 


SUBROUTINE  QPHESM  N,  NROHH,  NCOLH,  JTHCOL,  HESS,  X,  HX  ) 
INTEGER  N,  NROHH,  NCOLH,  JTHCOL 
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80  DOUBLE  PRECISION  HESS! NRONH, NCOLH),  HX(N),  X(N) 

C 

c  - 

C  OWES  I  COMPUTES  THE  VECTOR  HX  «  (HES3)*X  POR  SOME  MATRIX  HESS 
C  THAT  DEFINES  THE  HESSIAN  OF  THE  REQUIRED  QP  PROBLEM. 

C 

C  IN  THIS  VERSION  OF  QPHESS  THE  HESSIAN  MATRIX  IS  IMPLICIT. 

C  THE  ARRAY  HESS  IS  NOT  ACCESSED.  THERE  IS  NO  SPECIAL  COOINB 
C  FOR  THE  CASE  JTHCOL  .ST.  0. 

C  - 

c 

81  DOUBLE  PRECISION  ONE,  TNO 

82  DATA  ONE/1. 00*8/ ,  TNO/2.00+0/ 

C 

83  HX1 1)  =  THO*X<  1 1 

84  HX(2)  =  TNO*XC  2 ) 

85  HX(  3)  =  THO*(X(  3)  ♦  X(41) 

86  HX14)  s  HX(3I 

87  HX(5)  *  TM0*X(5) 

88  HX161  «  -  THO»(X(6)  ♦  X(7)1 

89  HX(  7)  s  HX(6) 

90  RETURN 
C 

C  END  OF  QPHES1 

91  END 

92  SUBROUTINE  qPHESZI  N.  NRONH,  NCOLM,  JTHCOL,  HESS,  X,  HX  1 

93  INTEGER  N,  NRONH,  NCOLH,  JTHCOL 

94  DOUBLE  PRECISION  HESS! NRONH, NCOLH),  HXIN),  X(N) 

C 

c  - - - - - 

C  IN  THIS  VERSION  OF  QPHESS.  THE  MATRIX  H  IS  STORED  IN  HESS  AS 
C  A  FULL  TNO-OIMENSIONAL  ARRAY. 

C  COPYVC  AND  ZEROVC  ARE  UTILITY  ROUTINES  USED  BY  QPSOL. 

C  - 

C 


95 

INTEGER 

I,  J 

96 

C 

DOUBLE  PRECISION 

XJ 

97 

C 

IF  (JTHCOL  .EQ.  0) 

GO  TO  100 

C  SPECIAL  CASE  —  EXTRACT  ONE  COLUFN  OF  H. 

C 

98  CALL  COPYVC1  N,  HES3( 1 .JTHCOL ) ,  N,  1,  HX,  N,  1  ) 

99  RETURN 
C 

C  NORMAL  CASE. 

C 

100  100  CALL  ZEROVCC  N,  HX,  N,  1  ) 

101  DO  200  J  »  1,  N 

102  XJ  *  X(J) 

103  00  150  I  «  1,  N 

104  HX(I)  *  HX(I)  ♦  HESS(I,J)*XJ 

105  150  CONTINUE 

106  200  CONTINUE 

107  RETURN 
C 

C  END  OF  QPHESE 
END 


108 
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109  SUBROUTINE  RPHES3(  N.  NROMG  NCOLM,  JTHCOL,  HESS.  X.  NX  ) 

110  INTEGER  N.  NROHH,  NCOLH.  JTHCOL 

111  DOUBLE  PRECISION  HESS! NR0M1  .NCOLH I >  HX(N).  XIN) 

C 

c  - 

C  IN  THIS  VERSION  OF  QPHESS,  THE  SYMMETRIC  PART  OF  H  IS  STORED  IN 
C  THE  LONER  HALF  OF  THE  TWO-DIMENSIONAL  ARRAY  HESS.  I.E.»  IN  THE 


c 

c 

ELEMENTS  HESS(I.J),  I  .6E.  J. 

112 

c 

INTEGER  I,  J,  JP1 ,  LRONH,  NM1 »  NUN 

113 

p 

DOUBLE  PRECISION  S.  XJ 

114 

p 

IF  (JTHCOL  .£«.  01  GO  TO  100 

V 

c 

r 

SPECIAL  CASE  —  EXTRACT  ONE  COLUMN  OF  H. 

115 

LRONH  i  NRONH<M  JTHCOL  -  1 1  ♦  1 

116 

CALL  COPYVCC  JTHCOL.  HESSC JTHCOL , 1 1 ,  LRONH, 

NROHH,  NX,  JTHCOL 

117 

NUH  s  N  -  JTHCOL 

118 

JP1  =  JTHCOL  ♦  1 

119 

IF  (HUM  .ST.  01 

•  CALL  COPYVCC  NUH,  HESS! JP1 .JTHCOL),  NUH, 

1,  HXIJPU, 

NUH, 

120 

p 

RETURN 

c 

p 

NORMAL  CASE. 

121 

100  DO  200  I  *  1,  N 

122 

S  =  0.0D«0 

123 

00  150  J  =  I,  N 

124 

S  =  S  ♦  HE33(J,I)«X(JI 

125 

150  CONTINUE 

126 

HX(I)  s  s 

127 

200  CONTINUE 

126 

p 

IF  (N  .LE.  1)  RETURN 

129 

NM1  *  N  -  1 

130 

DO  400  J  =  1 ,  NM1 

131 

XJ  =  X(J» 

132 

JP1  s  J  ♦  1 

133 

DO  350  I  s  JP1 ,  N 

134 

HX(I)  *  HX1 1 )  *  HESS(I,J)*XJ 

135 

350  CONTINUE 

136 

400  CONTINUE 

137 

p 

RETURN 

c 

END  OF  RPHES3 

138 

END 

139 

SUBROUTINE  RPHES41  N,  NROHH,  NCOLH,  JTHCOL, 

HESS,  X,  HX 

) 

140 

INTEGER  N,  NROHH,  NCOLH,  JTHCOL 

141 

DOUBLE  PRECISION  HESS! NROHH ) ,  HX(N),  X(N) 

C 

C  - 

C  IN  THIS  VERSION  OF  QPHESS,  THE  SYW1ETRIC  PART  OF  H  IS  STORED  IN 
C  THE  ONE-DIMENSIONAL  ARRAY  HESS.  NOTE  THAT  NROHH  IS  USED  TO  DEFINE 
C  THE  LENGTH  OF  HESS,  AND  MUST  BE  AT  LEAST  N*(N  ♦  1 )/2.  THE 
C  PARAMETER  NCOLH  IS  NOT  USED  HERE,  BUT  IT  MUST  BE  SET  TO  1  FOR 
C  THE  CALL  TO  RPSOL. 
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C 

142 

INTEGER  I*  INC.  J,  JPI,  L»  NKt.  NUN 

143 

DOUBLE  PRECISION  S.  XJ 

C 

144 

IF  (JTHCOL  .EG.  0)  60  TO  100 

t 

c 

SPECIAL  CASE  —  EXTRACT  ONE  COLIM4  OF  M. 

c 

145 

L  =  JTHCOL 

146 

INC  =  N 

147 

DO  50  I  «  1 .  JTHCOL 

148 

HX1I)  =  HESS1LI 

149 

INC  =  INC  -  1 

150 

L  =  L  ♦  INC 

151 

50  CONTINUE 

c 

152 

L  *  l  -  INC  ♦  1 

153 

NUN  a  N  -  JTHCOL 

154 

JPI  *  JTHCOL  ♦  1 

155 

IF  (HUM  .6T.  0) 

•  CALL  COPYVCt  NUN.  HESS(L).  NUN.  1.  HX(JPI),  NUN 

.  1  ) 

156 

RETURN 

c 

c 

NORMAL  CASE. 

c 

157 

100  L  a  0 

158 

DO  200  I  «  1.  N 

159 

S  *  0.00*0 

160 

DO  150  J  a  J,  N 

161 

L  a  t  ♦  1 

162 

S  a  S  ♦  HESS(L)«XU) 

163 

150  CONTINUE 

164 

HXIII  *  S 

165 

200  CONTINUE 

166 

f* 

IF  (N  .LE.  1 1  RETURN 

167 

Ir 

L  a  o 

168 

NMt  a  N  -  1 

169 

DO  400  J  a  | ,  NMt 

170 

XJ  a  X1J) 

171 

L  a  L  ♦  1 

172 

JPI  =  J  ♦  1 

173 

DO  350  I  a  Jpl,  N 

174 

L  a  L  ♦  1 

175 

HX(I)  a  HX(I)  *  HESS!  L)»XJ 

176 

350  CONTINUE 

177 

400  CONTINUE 

178 

RETURN 

1* 

C 

END  OF  RPHES4 

179 

END 

180 

SUBROUTINE  QPHE351  N.  NROMH.  NCOLH.  JTHCOL.  HESS, 

X,  HX  > 

181 

INTEGER  N,  NROMH,  NCOLH,  JTHCOL 

182 

DOUBLE  PRECISION  HESS1  NROMH, NCOLH),  KX(N),  XtN) 

C 

L 

c 

IN  THIS  VERSION  OF  GPHESS,  THE  CHOLE3KY  FACTOR  OF  N 

IS  STORED  IN 

c 

THE  LOWER  HALF  OF  THE  TWO-DIMENSIONAL  ARRAY  HESS.  IN  OTHER  WORDS; 

r  H  *  L  *  H TRANSPOSE).  WHERE  l  IS  A  LOWER  TRIANGULAR  MATRIX  'TOR ED 


v 
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a 


183 

184 

185 

186 


187 

188 

189 

190 

191 

192 


’  4 


C 

c 

c 


IN  HESSII.J),  I  .GE.  J. 


INTEGER 

INTEGER 

DOUBLE  PRECISION 


I.  IBACK,  J,  MAX,  LRONH,  NUN 

NINO 

S 


c 

c 

c 


IF  (JTHCOL  .EG.  0)  60  TO  100 

SPECIAL  CASE  --  HE  NEED  HX  =  L  «  ( JTH  RON  OF  L>. 

NUN  *  N  -  JTHCOL  ♦  1 

CALL  ZEROVC1  NUN,  HX( JTHCOL )>  NUN,  1  ) 

NUN  s  JTHCOL 

LROMH  s  NRONH* l NUM  -  11  ♦  1 

CALL  COPYVCt  NUN,  HESS! JTHCOL, 1 ),  LROMH,  NROW,  HX,  NUN,  t  > 
GO  TO  300 

NORMAL  CASE. 


193 

100  DO  200  I  S  1,  N 

* 

194 

S  =  0.00+0 

195 

DO  150  J  =  I,  N 

196 

S  =  S  ♦  HESS! J,I)*X(J) 

197 

150  CONTINUE 

198 

HXII)  =  S 

199 

200  CONTINUE 

r* 

200 

C 

NUM  =  N 

ft 

;rr 

C 

. 

C 

COMPUTE  HX  s  l  »  HX. 

,  - 

*J r 

c 

.  • 

201 

300  IBACK  =  N 

-  *■ 

•y  - 

202 

DO  400  I  s  1,  N 

• 

203 

S  =  0.00+0 

204 

JMAX  s  MINOI  NUM,  IBACK  1 

■ 

205 

00  350  J  s  1 »  JMAX 

206 

S  =  S  ♦  HESSI IBACK , J )*HXt  J ) 

*  v 

207 

350  CONTINUE 

208 

HXI  IBACK)  =  S 

209 

IBACK  *  IBACK  -  1 

210 

400  CONTINUE 

211 

RETURN 

w — V 

w 

c 

END  OF  GPHES5 

-  ; 

4".  * 

212 

END 

S'1 

'"0 

213 

SUBROUTINE  QPHES6I  N,  NRONH,  NCOLH,  JTHCOL,  HESS,  X,  HX  1 

*«  •«* 
w. « 

‘  •  \ 

214 

INTEGER  N,  NRONH,  NCOLH,  JTHCOL 

,  - 

•■>,  4 

215 

DOUBLE  PRECISION  HESSI NRONH),  HXIN),  XI N) 

•  •  * 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


IN  THIS  VERSION  OF  GPHESS,  THE  CHOLE3KY  FACTOR  OF  H  IS  STORED  IN 
THE  ONE-DIMENSIONAL  ARRAY  HESS.  IN  OTHER  WORDS, 

H  =  L  »  LI  TRANSPOSE),  WERE  L  IS  A  LONER  TRIANGULAR  MATRIX  STORED 
COMPACTLY  BY  COLUMNS  IN  HESS.  NOTE  THAT  NRONH  IS  USEO  TO  DEFINE 
THE  LENGTH  OF  HESS,  AND  MUST  BE  AT  LEAST  WIN  ♦  1)/2.  THE 
PARAMETER  NCOLH  IS  NOT  USED  HERE,  BUT  IT  SHOULD  BE  SET  TO  1  FOR 
THE  CALL  TO  QP30L. 


IVJ 

<•  •  m  • . 
*>  r  V 

v 


,  n 


A.vV\A.S.S  .  ..N  t .....  . 


•  W.V-VV. ; 
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CVEC. 


-0.02 

-0.20 

-0.20  -0.20 

-0.20 

0.04 

0.04 

ROUS  OF  A. 

1.00 

1.00 

1.00  1.00 

1.00 

1.00 

1.00 

0.15 

0.04 

0.02  0.04 

0.02 

0.01 

0.03 

0.03 

0.05 

0.08  0.02 

0.06 

0.01 

0.00 

0.02 

0.04 

0.01  0.02 

0.02 

0.00 

e.oo 

0.02 

0.03 

0.00  0.00 

0.01 

0.00 

0.00 

0.70 

0.75 

0.80  0.75 

0.80 

0.97 

0.00 

0.02 

0.06 

0.08  0.12 

0.02 

0.01 

0.97 

LOUER  BOUNDS. 

-0.100-01  -0. 

10D  00  -0.100-01  -0.400-01 

-0.100  00 

-0.10D-01 

-0.10D-0I 

-0.130  00  -0. 

100  13  - 

■0.100  13  -0.1  OD  13 

-0.10D  13 

-0.990-01 

-0.30D-02 

UPPER  BOUNDS. 

0.100-01  0. 

150  00 

0.300-01  0.200-01 

0.50D-01 

0.100  13 

0.10D  13 

-0.130  00  -0. 

490-02  - 

-0.64D-02  -0.370-02 

-0. 12D-02 

0.100  13 

0.20D-02 

INITIAL  X. 

-0.01 

-0.03 

0.00  -0.01 

-0.10 

0.02 

0.01 

M0RK3PACE  PROVIDED  IS 

IH(  501,  HI 

200). 

TO  SOLVE  PROBLEM  HE  NEED  IM(  141,  HI 

161). 

ITN  JOEL  JAOD  STEP 

000  0.00D-01 

1  W  tJl  4.120-02 

2  12U  41  4.240-02 


COND  T  NUHINF 
1.030  02  3 
1.560  02  1 
5.300  01  0 


SUMINF 

1.0300000-01 

3.0000000-02 

0.0000000-01 


EXIT  LP  PHASE.  INFORM  *  0  ITER  =  2 


ITN 

JOEL 

JAOD 

STEP 

NHES3 

OBJECTIVE 

NCOLZ  NORM  6FREE 

NORM  ZT6 

COND  T  COND  ZHZ 

HESS  MOO 

0 

0 

0 

0.000-01 

1 

4.58000-02 

0 

2.410-01 

0.000-01 

5.30  01 

1.00 

00 

0.000-01 

0 

5L 

0 

O.OOD-OI 

2 

4.58000-02 

1 

4.670-01 

2. 160-01 

6.0D  01 

1.00 

00 

0.000-01 

1 

0 

14L 

1.330-01 

3 

4.16160-02 

0 

4.440-01 

O.OOD-OI 

6.00  01 

1.00 

00 

0.000-01 

1 

11U 

0 

0.000-01 

4 

4.16160-02 

1 

4.440-01 

9.460-02 

1.30  01 

1.00 

00 

0.000-01 

2 

0 

0 

1.000  00 

5 

3.93620-02 

t 

4.330-01 

1.390-17 

1.30  01 

1 .00 

00 

0.000-01 

2 

3L 

0 

0.000-01 

6 

3.93620-02 

2 

5.260-01 

9.200-02 

1.50  01 

1.30 

00 

0.000-01 

3 

0 

10U 

4.150-01 

7 

3.75890-02 

1 

5.180-01 

1.190-02 

5.70  01 

1.00 

00 

0.000-01 

4 

0 

0 

1 .000  00 

8 

3.75540-02 

1 

5.180-01 

3.47D-18 

5.70  01 

1.00 

00 

0.000-01 

4 

4L 

0 

0.000-01 

9 

3.75540-02 

2 

5.770-01 

5.010-02 

5.30  01 

1 .20 

00 

0.000-01 

5 

0 

0 

1.000  00 

10 

3.70320-02 

2 

5.570-01 

8.590-18 

5.30  01 

1 .20 

00 

0.000-01 

EXIT  OP  PHASE. 

INFORM  =  0  ITER  = 

5 

VARIABLE 

STATE 

VALUE 

LONER  BOUND 

UPPER  BOUND 

LA6R  MULTIPLIER 

RESIDUAL 

VARBL 

1 

LL 

-0.10000000-01 

-0.10000000-01 

0.10000000-01 

0.4700306 

0.0000 

VARBL 

2 

FR 

-0.69864650-01 

-0.1000000 

0.1500000 

0.0000000 

0.30140-01 

VARBL 

3 

FR 

0.18259150-01 

-0. 10000000-01 

0.30000000-01 

0.0000000 

0.11740-01 

VARBL 

4 

FR 

-0.24260810-01 

-0.40000000-01 

0.20000000-01 

0.0000000 

0.15740-01 

VARBL 

5 

FR 

-0.62005640-01 

-0.1000000 

0.50000000-01 

0.0000000 

0.37990-01 

VARBL 

6 

FR 

0.13805440-01 

-0.10000000-01 

NONE 

o.r"ooooo 

0.23810-01 

VARBL 

7 

FR 

0.40664960-02 

-0.10000000-01 

NONE 

0.0000000 

0.1407D-0I 
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LINEAR 

CONSTR 

STATE 

VALUE 

LONER  BOUIO 

UPPER  BOUND  LAGR  MULTIPLIER 

RESIDUAL 

LNCON 

1 

EQ 

-0.1300000 

-0.1300000 

-0.1300000 

-1.908183 

0.41630-16 

LNCON 

2 

FR 

-0.58798980-02 

NONE 

-0.49000000-02 

0.0000000 

0.97990-03 

LNCON 

3 

UL 

-0. 64000000-02 

NONE 

-0.64000000-02 

-0.3143604 

0.86740-18 

LNCON 

4 

FR 

-0.45373230-02 

NONE 

-0.37000000-02 

0.0000000 

0.83730-03 

LNCON 

5 

FR 

-0.2915996D-02 

NONE 

-0.12000000-02 

0.0000000 

0.17160-02 

LNCON 

6 

LL 

-0. 99200000-01 

-0.99200000-01 

NONE 

1.954501 

0.55510-16 

LNCON 

7 

LL 

-0.30000000-02 

-0.30000000-02 

0.20000000-02 

1.971586 

0.27110-18 

EXIT  QPSOL  -  OPTIMAL  QP  SOLUTION. 

FINAL  qP  OBJECTIVE  VALUE  -  0.37031650-01 


A  RUN  OF  THE  SAME  EXAMPLE  WITH  A  HARM  START.... 

WORKSPACE  PROVIDED  IS  IM(  501.  H(  200). 

TO  SOLVE  PROBLEM  ME  NEED  IH1  14),  H(  161). 

EXIT  LP  PHASE.  INFORM  *  0  ITER  =  0 


ITN  JOEL  JADO  STEP  NHESS  OBJECTIVE  NCOLZ  NORM  6FREE  NORM  ZTG  COMO  T  COND  2HZ  HESS  MOD 

000  0.00D-01  3  3.70320-02  2  5.57D-01  8.650-16  3.50  01  1.30  00  0.000-01 

EXIT  QP  PHASE.  INFORM  =  0  ITER  =  0 

EXIT  QPSOL  -  OPTIMAL  QP  SOLUTION. 


FINAL  QP  OBJECTIVE  VALUE  a  0.37031650-01 
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so.  abstract  re— ibnm  —  RWRt  Mas  u  b»»« nay  mt  i— bp  m— *  — This  report  forms  the 
user's  guide  for  Version  3.2  of  QPSOL,  a  set  of  Fortran  subroutines  designed  to 
locate  the  minimum  value  of  an  arbitrary  quadratic  function  subject  to  linear 
constraints  and  simple  upper  and  lower  bounds.  If  the  quadratic  function  is 
convex,  a  global  minimum  is  found;  otherwise,  a  local  minimum  is  found.  The 
method  used  is  most  efficient  sfcen  many  constraints  or  bounds  are  active  at  the 
solution.  QPSOL  treats  the  Hessian  and  general  constraints  as  dense  matrices, 
and  hence  is  not  intended  for  large  sparse  problems.  This  document  replaces  the 
previous  user's  guide  of  Duly  1983. 
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